Tom LaSalle
This document contains all the code necessary to generate the plots for Figure 1 and related supplementary figures (S1-S6). Plots are subsequently edited in Adobe Illustrator to produce the final figures.
Load the necessary libraries:
library(knitr)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)
library(plyr)
library(dplyr)
library(DESeq2)
library(openxlsx)
library(corrplot)
library(ggmosaic)
library(cowplot)
library(Seurat)
library(heatmap3)
library(ggpubr)
library(umap)
library(fgsea)
Import the metadata and keep only data for which neutrophil enrichment was performed:
prefix <- "~/Downloads/Github/"
metadata_bypatient <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 2)
qc_data <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 9)
metadata_bypatient <- metadata_bypatient[which(metadata_bypatient$Public.ID %in% qc_data$Public.ID),]
genomic_signatures <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 10)
metadata_long <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 4)
Due to an institutional IRB-approved waiver of informed consent, all clinical data are reported in quintiles, and not all clinical characteristics are publicly available. Correlations were calculated between the following clinical variables: Age, Sex, BMI, Heart condition, Lung condition, Kidney condition, Diabetes, Hypertension, Immunocompromised, Symptom duration, Respiratory symptoms, Fever symptoms, GI Symptoms, Acuity (D0, D3, D7, Max within D7-D28, D28, Max Overall), ANC (D0, D3, D7), ALC (D0, D3, D7), AMC (D0, D3, D7), Creatinine (D0, D3, D7), CRP (D0, D3, D7), D-dimer (D0, D3, D7), LDH (D0, D3, D7), Troponin detected within 72 hours, Intubation status, Race (White, Black, Asian, Other), Hispanic ethnicity, CXR infiltrates, and Smoking history. Multiple hypothesis correction was performed to determine significant correlations, and then those parameters with significant correlations were plotted in a correlation heatmap. Though it is not possible to show all of the data, the code to generate the figure is below and is written as if the additional columns were present in the metadata file:
# # Reverse factor levels of Acuity such that positive correlations with Acuity indicate correlations with worse disease severity
# metadata_bypatient$Acuity.0 <- mapvalues(metadata_bypatient$Acuity.0, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.3 <- mapvalues(metadata_bypatient$Acuity.3, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.7 <- mapvalues(metadata_bypatient$Acuity.7, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.28 <- mapvalues(metadata_bypatient$Acuity.28, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.7.to.28 <- mapvalues(metadata_bypatient$Acuity.7.to.28, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.max <- mapvalues(metadata_bypatient$Acuity.max, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
#
# rownames(metadata_bypatient) <- metadata_bypatient$Public.ID
#
# metadata_bypatient <- metadata_bypatient[,-which(colnames(metadata_bypatient) %in% c("Study_ID","Public.ID","D0_draw","D3_draw","D7_draw","DE_draw","D0_Neu_Sample","D3_Neu_Sample","D7_Neu_Sample","DE_Neu_Sample"))]
# metadata_bypatient <- metadata_bypatient[metadata_bypatient$COVID == "1",]
# metadata_bypatient <- metadata_bypatient[,-which(colnames(metadata_bypatient) %in% c("COVID"))]
#
# cormatrix <- matrix(0L, nrow = ncol(metadata_bypatient), ncol = ncol(metadata_bypatient))
# pmatrix <- matrix(0L, nrow = ncol(metadata_bypatient), ncol = ncol(metadata_bypatient))
# rownames(cormatrix) <- rownames(pmatrix) <- colnames(cormatrix) <- colnames(pmatrix) <- colnames(metadata_bypatient)
#
# for (i in 1:nrow(cormatrix)){
# for (j in 1:nrow(cormatrix)){
# stats <- cor.test(x = as.numeric(metadata_bypatient[,i]), y = as.numeric(metadata_bypatient[,j]), use = "pairwise.complete.obs", method = "kendall")
# cormatrix[i,j] <- stats$estimate
# pmatrix[i,j] <- stats$p.value
# }
# }
cormatrix <- as.matrix(read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 5, rowNames = TRUE))
pmatrix <- as.matrix(read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 6, rowNames = TRUE))
ptable <- matrix(0L, nrow = nrow(pmatrix)^2, ncol = 3)
colnames(ptable) <- c("Row","Column","pvalue")
ptable[,1] <- rep(colnames(pmatrix), each = nrow(pmatrix))
ptable[,2] <- rep(colnames(pmatrix), nrow(pmatrix))
for (i in 1:nrow(ptable)){
ptable[i,3] <- pmatrix[rownames(pmatrix) %in% ptable[i,1],colnames(pmatrix) %in% ptable[i,2]]
}
ptable <- as.data.frame(ptable[complete.cases(ptable[,3]),])
ptable$fdr <- p.adjust(ptable[,3], method = "fdr")
ptable <- ptable[ptable[,4] < 0.05,]
ptable <- ptable[ptable[,1] %in% c("abs_neut_0_cat","abs_neut_3_cat","abs_neut_7_cat") | ptable[,2] %in% c("abs_neut_0_cat","abs_neut_3_cat","abs_neut_7_cat"),]
cormatrix <- cormatrix[rownames(cormatrix) %in% unique(c(ptable[,1],ptable[,2])),colnames(cormatrix) %in% unique(c(ptable[,1],ptable[,2]))]
pmatrix <- pmatrix[rownames(pmatrix) %in% unique(c(ptable[,1],ptable[,2])),colnames(pmatrix) %in% unique(c(ptable[,1],ptable[,2]))]
rownames(cormatrix) <- colnames(cormatrix) <- rownames(pmatrix) <- colnames(pmatrix) <- mapvalues(rownames(cormatrix), from = c("BMI.cat","abs_mono_0_cat","abs_mono_3_cat","abs_mono_7_cat","Age.cat","creat_0_cat","creat_3_cat","creat_7_cat","ldh_0_cat","ldh_3_cat","ldh_7_cat","crp_0_cat","crp_3_cat","crp_7_cat","Acuity.28","Acuity.7","Acuity.3","Acuity.0","Acuity.max","abs_neut_0_cat","abs_neut_3_cat","abs_neut_7_cat","ddimer_0_cat","ddimer_3_cat","ddimer_7_cat"), to = c("BMI","AMC, D0", "AMC, D3","AMC, D7","Age","Creatinine, D0","Creatinine, D3","Creatinine, D7","LDH, D0","LDH, D3","LDH, D7","CRP, D0","CRP, D3","CRP, D7","Acuity, D28","Acuity, D7","Acuity, D3","Acuity, D0","Acuity, Max","ANC, D0","ANC, D3","ANC, D7","D-dimer, D0","D-dimer, D3","D-dimer, D7"))
Figure 1B:
corrplot(cormatrix, method = "square", type = "lower", order = "hclust", hclust.method = "ward.D", tl.col="black", col=colorRampPalette(c("blue","white","red"))(200), p.mat = pmatrix, sig.level = 0.05, insig = "blank")
Next we highlight the correlations between ANC on each individual day and Acuity Max, using Kendall’s tau correlation:
metadata_long <- metadata_long[which(metadata_long$Public.ID %in% qc_data$Public.ID),]
metadata_temp <- metadata_long[complete.cases(metadata_long$ANC.matched) & metadata_long$COVID == 1,]
metadata_temp$Acuity.max <- factor(metadata_temp$Acuity.max, levels = c("5","4","3","2","1"))
corvals0 <- cor.test(x = as.numeric(metadata_temp[metadata_temp$Day == "D0",]$Acuity.max), y = as.numeric(metadata_temp[metadata_temp$Day == "D0",]$ANC.matched), use = "pairwise.complete.obs", method = "kendall")
corvals3 <- cor.test(x = as.numeric(metadata_temp[metadata_temp$Day == "D3",]$Acuity.max), y = as.numeric(metadata_temp[metadata_temp$Day == "D3",]$ANC.matched), use = "pairwise.complete.obs", method = "kendall")
corvals7 <- cor.test(x = as.numeric(metadata_temp[metadata_temp$Day == "D7",]$Acuity.max), y = as.numeric(metadata_temp[metadata_temp$Day == "D7",]$ANC.matched), use = "pairwise.complete.obs", method = "kendall")
metadata_temp$ANC.matched <- factor(metadata_temp$ANC.matched, levels = c("5","4","3","2","1"))
my.cols <- brewer.pal(9, "RdBu")
p1 <- ggplot(data = metadata_temp[metadata_temp$Day == "D0",]) + geom_mosaic(aes(x = product(ANC.matched, Acuity.max), fill = ANC.matched), offset = 0.005) + coord_fixed(ratio = 1) + theme_bw() + scale_fill_manual(values = (my.cols[(c(1,2,3,7,9))])) + ggtitle("Day 0") + annotate(geom="text", x=0.5, y=.8, label=paste0("τ =",round(cormatrix[rownames(cormatrix) == "ANC, D0",colnames(cormatrix) == "Acuity, Max"],3),"\npadj =",round(pmatrix[rownames(pmatrix) == "ANC, D0",colnames(pmatrix) == "Acuity, Max"],12)), color="black") + theme(legend.key.size=unit(4,'mm'), legend.text=element_text(size=5), legend.title=element_text(size=5), axis.title = element_text(size = 7))
p2 <- ggplot(data = metadata_temp[metadata_temp$Day == "D3",]) + geom_mosaic(aes(x = product(ANC.matched, Acuity.max), fill = ANC.matched), offset = 0.005) + coord_fixed(ratio = 1) + theme_bw() + scale_fill_manual(values = (my.cols[(c(1,2,3,7,9))])) + ggtitle("Day 3") + annotate(geom="text", x=0.5, y=.8, label=paste0("τ =",round(cormatrix[rownames(cormatrix) == "ANC, D3",colnames(cormatrix) == "Acuity, Max"],3),"\npadj =",round(pmatrix[rownames(pmatrix) == "ANC, D3",colnames(pmatrix) == "Acuity, Max"],11)), color="black") + theme(legend.key.size=unit(4,'mm'), legend.text=element_text(size=5), legend.title=element_text(size=5), axis.title = element_text(size = 7))
p3 <- ggplot(data = metadata_temp[metadata_temp$Day == "D7",]) + geom_mosaic(aes(x = product(ANC.matched, Acuity.max), fill = ANC.matched), offset = 0.005) + coord_fixed(ratio = 1) + theme_bw() + scale_fill_manual(values = (my.cols[(c(1,2,3,7,9))])) + ggtitle("Day 7") + annotate(geom="text", x=0.5, y=.8, label=paste0("τ =",round(cormatrix[rownames(cormatrix) == "ANC, D7",colnames(cormatrix) == "Acuity, Max"],3),"\npadj =",round(pmatrix[rownames(pmatrix) == "ANC, D7",colnames(pmatrix) == "Acuity, Max"],11)), color="black") + theme(legend.key.size=unit(4,'mm'), legend.text=element_text(size=5), legend.title=element_text(size=5), axis.title = element_text(size = 7))
Figure 1C:
plot_grid(p1,p2,p3,ncol=3)
Now we begin to look at our transcriptomic data. We import the count and TPM data from RSEM, as well as the Ensembl ID to gene symbol conversions, and log-transform the TPM values:
Count <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 7, rowNames = TRUE)
TPM <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 8, rowNames = TRUE)
genepc <- read.delim(paste0(prefix,"Ensembl_to_Symbol.txt"))
logTPM <- log2(TPM + 1)
metadata_long <- merge(metadata_long, qc_data)
rownames(metadata_long) <- metadata_long$Public.Sample.ID
To get an overall sense of the data we start with a PCA of the protein-coding genes:
proteincoding <- genepc$Gene.stable.ID[genepc$Gene.type == "protein_coding"]
logPCG <- logTPM[which(rownames(logTPM) %in% proteincoding),]
logPCG <- logPCG[which(rowSums(logPCG) > 0),]
pcaResults <- prcomp(t(logPCG), center = T, scale. = T)
PCs <- data.frame(pcaResults$x[,1:2])
metadata_long$PC1 <- PCs$PC1
metadata_long$PC2 <- PCs$PC2
The main contributions to PC1 are Exonic rate and Median 3’ bias:
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(0,70))
p1 <- ggplot(data = metadata_long, aes(x = -1*PC1, y = PC2, color = Exonic.Rate*100)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank()) + ggtitle("PCA") + sc + coord_fixed(ratio = 2.2)
p1$labels$colour <- "Exonic Rate"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlGn")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_long$Median.3..bias),1))
p2 <- ggplot(data = metadata_long, aes(x = -1*PC1, y = PC2, color = Median.3..bias)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank()) + ggtitle("PCA") + sc + coord_fixed(ratio = 2.2)
p2$labels$colour <- "Median 3' Bias"
Figure S1B:
plot_grid(p1,p2)
We apply quality control filtration to the data and re-examine the PCA:
#Calculate mitochondrial percentages
mitorows <- grepl("^MT-",genepc$Gene.name)
mitoids <- genepc$Gene.stable.ID[mitorows]
mitocounts <- Count[rownames(Count) %in% mitoids,]
mitosums <- colSums(mitocounts)
totalcounts <- colSums(Count)
percent.mt <- mitosums/totalcounts*100
metadata_long$percent.mt <- percent.mt
#Perform quality control filtration based on RNASeQC parameters
metadata_filtered <- metadata_long[metadata_long$percent.mt < 20 & metadata_long$Genes.Detected > 10000 & metadata_long$Median.Exon.CV < 1 & metadata_long$Exon.CV.MAD < 0.75 & metadata_long$Exonic.Rate*100 > 25 & metadata_long$Median.3..bias < 0.9,]
logPCG_filtered <- logPCG[,colnames(logPCG) %in% rownames(metadata_filtered)]
logTPM_filtered <- logTPM[,colnames(logTPM) %in% rownames(metadata_filtered)]
TPM_filtered <- TPM[,colnames(TPM) %in% rownames(metadata_filtered)]
Count_filtered <- Count[,colnames(Count) %in% rownames(metadata_filtered)]
tf <- rowSums(TPM_filtered > 0.1) > ncol(TPM_filtered)*.2
TPM_filtered <- TPM_filtered[tf,]
Count_filtered <- Count_filtered[tf,]
logTPM_filtered <- logTPM_filtered[tf,]
tf <- rowSums(Count_filtered >= 6) > ncol(Count_filtered)*.2
TPM_filtered <- TPM_filtered[tf,]
Count_filtered <- Count_filtered[tf,]
logTPM_filtered <- logTPM_filtered[tf,]
#Repeat PCA
logPCG_filtered <- logPCG_filtered[which(rowSums(logPCG_filtered) > 0),]
pcaResults <- prcomp(t(logPCG_filtered), center = T, scale. = T)
PCs <- data.frame(pcaResults$x[,1:2])
metadata_filtered$PC1 <- PCs$PC1
metadata_filtered$PC2 <- PCs$PC2
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(0,70))
p1 <- ggplot(data = metadata_filtered, aes(x = -1*PC1, y = -1*PC2, color = Exonic.Rate*100)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank()) + ggtitle("PCA") + sc + coord_fixed(ratio = 1.3)
p1$labels$colour <- "Exonic Rate"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlGn")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered$Median.3..bias),1))
p2 <- ggplot(data = metadata_filtered, aes(x = -1*PC1, y = -1*PC2, color = Median.3..bias)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank()) + ggtitle("PCA") + sc + coord_fixed(ratio = 1.3)
p2$labels$colour <- "Median 3' Bias"
Figure S1C:
plot_grid(p1,p2)
We check that quality control filtration has not introduced any biases in terms of patient demographics.
metadata_long$QC <- rownames(metadata_long) %in% rownames(metadata_filtered)
metadata_long$QC <- mapvalues(metadata_long$QC, from = c("FALSE","TRUE"), to = c("FAIL","PASS"))
histmaker <- function(parameter){
metadata_temp <- metadata_long[,colnames(metadata_long) %in% c(parameter,"QC")]
metadata_temp[,1] <- factor(metadata_temp[,1])
histtable <- as.data.frame(matrix(0L, nrow = nlevels(metadata_temp[,1])*2, ncol = 3))
colnames(histtable) <- c("Param","QC","Freq")
histtable[,1] <- rep(levels(metadata_temp[,1]), each = 2)
histtable[,2] <- rep(c("FAIL","PASS"),nlevels(metadata_temp[,1]))
for (i in 1:nrow(histtable)){
histtable$Freq[i] <- sum(metadata_temp[,1] == histtable$Param[i] & metadata_temp$QC == histtable$QC[i])
}
return(histtable)
}
parameter <- "COVID"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p1 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .0033) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
parameter <- "Acuity.max"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p2 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .018) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
parameter <- "severity.max"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p3 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
parameter <- "Age"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p4 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .024) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
#parameter <- "sex"
#fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
#histtable <- histmaker(parameter)
#p5 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
p5 <- ""
#parameter <- "race"
#fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
#histtable <- histmaker(parameter)
#p6 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
p6 <- ""
#parameter <- "ethnicity"
#fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
#histtable <- histmaker(parameter)
#p7 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
p7 <- ""
parameter <- "BMI"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p8 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .024) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
Figure S1E:
p1
p2
p3
p4
p8
We do observe that more samples are discarded from later time points:
parameter <- "Day"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p1 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .014) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: p = ",round(fishertest$p.value, digits = 8)), colour = "black")
Figure S1F:
p1
This may be partially attributed to longer sample processing times.
Figure S1G:
wilcoxtest <- wilcox.test(Processing.Time ~ QC, metadata_long)
ggplot(metadata_long, aes(x = factor(QC), y = as.numeric(Processing.Time), fill = factor(QC))) + geom_boxplot(outlier.shape = NA) + geom_jitter() + theme_bw() + ylab("Processing Time (h)") + scale_fill_manual(values = c("darkred","forestgreen")) + coord_fixed(ratio = .08) + annotate("text", x=1, y=23, label= paste0("Wilcoxon: p = ",round(wilcoxtest$p.value, digits = 5)), colour = "black")
However, there are no significant differences in processing time between samples that pass or fail within the same day.
Figure S1H:
ggplot(metadata_long, aes(x = factor(Day), y = Processing.Time, fill = QC)) + geom_boxplot(outlier.shape = NA) + geom_point(position = position_jitterdodge(jitter.height = 0, jitter.width = 0.2), size = 0.5, alpha = .5) + theme_bw() + scale_fill_manual(values = c("darkred","forestgreen")) + xlab("Day") + ylab("Processing Time (h)") + coord_fixed(ratio = 0.21)
Next we want to get an estimate of the sample purity since we performed bulk RNA-seq of enriched neutrophils. In order to do so, we need a reference single-cell RNA-seq dataset for deconvolution using CIBERSORTx. We will use the Bonn Cohort 2 from Schulte-Schrepping et al. We start by importing the data and collapsing the cluster labels into major hematopoeitic cell lineages.
seuratwb <- readRDS(paste0(prefix,"seurat_COVID19_freshWB-PBMC_cohort2_rhapsody_jonas_FG_2020-08-18.rds"))
barcodecluster <- as.character(seuratwb$RNA_snn_res.0.8)
temp <- seuratwb$RNA_snn_res.0.8
names(barcodecluster) <- names(temp)
barcodecluster <- cbind(rownames(barcodecluster),barcodecluster)
barcodecluster <- as.data.frame(barcodecluster)
barcodecluster$barcode <- as.character(rownames(barcodecluster))
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("0","1","4","8")] <- "Mature_Neutrophil"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("12","15")] <- "Immature_Neutrophil"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("3","10","14","18")] <- "Monocyte"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("2","5","6","13","16")] <- "T_NK"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("7","22")] <- "B"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("19")] <- "Plasmablast"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("9","11","17","20","21","23","24")] <- "Other"
seuratwb <- AddMetaData(seuratwb, metadata = subset(barcodecluster, select = c("barcode")), col.name = "Cell.Type")
Figure S2A:
DimPlot(seuratwb, reduction = "umap", group.by = "Cell.Type", label = FALSE) + scale_color_manual(values = c("grey60","tomato","forestgreen","skyblue","grey90","black","slateblue3")) + theme(axis.line.x = element_blank(), axis.line.y = element_blank(), axis.title.x = element_blank(), axis.ticks = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank())
For CIBERSORTx, we will generate a cell type signature matrix based on “pseudobulk” samples from a single cell type. We will add up the counts from all the cells of a given lineage for each patient and treat them as bulk samples.
pseudobulk <- matrix(0L, nrow = nrow(seuratwb@assays$RNA@counts), ncol = (length(levels(factor(seuratwb$donor))))*6)
rownames(pseudobulk) <- (seuratwb@assays$RNA@data@Dimnames[[1]])
idx <- (levels(factor(seuratwb$donor)))
pseudobulk <- as.data.frame(pseudobulk)
for (i in 0:(ncol(pseudobulk)/6 - 1)){
k <- i*6+1
colnames(pseudobulk)[k] <- paste(as.character(idx[i+1]),"_Mature_Neu",sep = "")
colnames(pseudobulk)[k+1] <- paste(as.character(idx[i+1]),"_Immature_Neu",sep = "")
colnames(pseudobulk)[k+2] <- paste(as.character(idx[i+1]),"_Mono",sep = "")
colnames(pseudobulk)[k+3] <- paste(as.character(idx[i+1]),"_TNK",sep = "")
colnames(pseudobulk)[k+4] <- paste(as.character(idx[i+1]),"_B",sep = "")
colnames(pseudobulk)[k+5] <- paste(as.character(idx[i+1]),"_Plasmablast",sep = "")
}
#Patient BN-31 does not have all the cell types so they are excluded
for (i in 0:17){
k <- i*6+1
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Mature_Neutrophil")
pseudobulk[,k] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Immature_Neutrophil")
pseudobulk[,k+1] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Monocyte")
pseudobulk[,k+2] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "T_NK")
pseudobulk[,k+3] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "B")
pseudobulk[,k+4] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Plasmablast")
pseudobulk[,k+5] <- rowSums(temp@assays$RNA@counts)
}
for (i in 19:(ncol(pseudobulk)/6 - 1)){
k <- i*6+1
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Mature_Neutrophil")
pseudobulk[,k] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Immature_Neutrophil")
pseudobulk[,k+1] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Monocyte")
pseudobulk[,k+2] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "T_NK")
pseudobulk[,k+3] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "B")
pseudobulk[,k+4] <- rowSums(temp@assays$RNA@counts)
temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Plasmablast")
pseudobulk[,k+5] <- rowSums(temp@assays$RNA@counts)
}
pseudobulk <- pseudobulk[rowSums(pseudobulk) > 0,]
pseudobulk <- pseudobulk[-grep("\\.",rownames(pseudobulk)),]
pseudobulk <- pseudobulk[-grep("^MT",rownames(pseudobulk)),]
Then we convert the Ensembl gene IDs to symbols for the neutrophil data and filter the pseudobulk matrix to include only the overlapping genes between the two datasets:
genepctemp <- genepc[genepc$Gene.stable.ID %in% rownames(Count_filtered),]
Count_temp <- Count
for (i in 1:nrow(Count_temp)){
if (length(which(genepctemp$Gene.stable.ID == rownames(Count_temp)[i])) == 1){
if (!(genepctemp$Gene.name[which(genepctemp$Gene.stable.ID == rownames(Count_temp)[i])] %in% rownames(Count_temp))){
rownames(Count_temp)[i] <- genepctemp$Gene.name[which(genepctemp$Gene.stable.ID == rownames(Count_temp)[i])]
}
}
}
pseudobulk_temp <- pseudobulk[rownames(pseudobulk) %in% rownames(Count_temp),]
The matrix that was used to generate the signature matrix is included in the Github.
Now we are able to run CIBERSORTx Generate Signature Matrix, setting limits of 50 to 100 genes per cell type and filtering for only hematopoietic genes. The CIBERSORTx cell type signature matrix:
Figure S2B:
matrix <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 13, rowNames = TRUE)
par(mar=c(1,1,1,1))
heatmap3(matrix, Colv = NA, cexRow = 0.2)
Using this cell type signature matrix, we generate estimated cell type fractions in CIBERSORTx using default parameters. An overview of the cell type proportions:
rownames(genomic_signatures) <- genomic_signatures$Public.Sample.ID
cibersort_temp <- genomic_signatures[,colnames(genomic_signatures) %in% c("Mature_Neutrophil","Immature_Neutrophil","Monocyte","T_NK","B","Plasmablast","Neutrophil_total")]
cibersort_temp <- cibersort_temp[rownames(cibersort_temp) %in% rownames(metadata_filtered),]
cibersort_sorted <- cibersort_temp[order(cibersort_temp$Neutrophil_total),]
cibersort_sorted <- cibersort_sorted[,-7]
data_percentage <- t(cibersort_sorted*100)
coul <- c("forestgreen","tomato","skyblue","slateblue3","gray","black")
Figure S2C:
barplot(data_percentage, col=coul , border=NA, xlab=NA, axisnames = FALSE)
Search for differences in CIBERSORTx estimated cell fractions across COVID status for all samples:
metadata_filtered <- merge(metadata_filtered, genomic_signatures)
rownames(metadata_filtered) <- metadata_filtered$Public.Sample.ID
metadata_filtered$COVID <- mapvalues(metadata_filtered$COVID, from = c(0,1), to = c("Negative","Positive"))
my.cols <- brewer.pal(3, "Set2")
wilcoxtest <- wilcox.test(as.numeric(Neutrophil_total)*100 ~ factor(COVID), metadata_filtered)
p1 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Neutrophil_total)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("Total")
wilcoxtest <- wilcox.test(as.numeric(Mature_Neutrophil)*100 ~ factor(COVID), metadata_filtered)
p2 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Mature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Mature")
wilcoxtest <- wilcox.test(as.numeric(Immature_Neutrophil)*100 ~ factor(COVID), metadata_filtered)
p3 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Immature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Immature")
wilcoxtest <- wilcox.test(as.numeric(Monocyte)*100 ~ factor(COVID), metadata_filtered)
p4 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Monocyte)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Monocyte")
wilcoxtest <- wilcox.test(as.numeric(T_NK)*100 ~ factor(COVID), metadata_filtered)
p5 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(T_NK)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 4)), colour = "black") + ggtitle("T/NK")
wilcoxtest <- wilcox.test(as.numeric(B)*100 ~ factor(COVID), metadata_filtered)
p6 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(B)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("B")
wilcoxtest <- wilcox.test(as.numeric(Plasmablast)*100 ~ factor(COVID), metadata_filtered)
p7 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Plasmablast)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 10)), colour = "black") + ggtitle("Plasmablast")
Figure S2D:
plot_grid(p1,p2,p3,p4,p5,p6,p7,ncol = 7)
Search for differences in CIBERSORTx estimated cell fractions across disease severity for all samples regardless of COVID status:
metadata_filtered <- merge(metadata_filtered, genomic_signatures)
metadata_temp <- metadata_filtered[metadata_filtered$severity.max %in% c("severe","non-severe"),]
my.cols <- brewer.pal(3, "RdBu")
wilcoxtest <- wilcox.test(as.numeric(Neutrophil_total)*100 ~ factor(severity.max), metadata_temp)
p1 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Neutrophil_total)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 18)), colour = "black") + ggtitle("Total")
wilcoxtest <- wilcox.test(as.numeric(Mature_Neutrophil)*100 ~ factor(severity.max), metadata_temp)
p2 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Mature_Neutrophil)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("Mature")
wilcoxtest <- wilcox.test(as.numeric(Immature_Neutrophil)*100 ~ factor(severity.max), metadata_temp)
p3 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Immature_Neutrophil)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Immature")
wilcoxtest <- wilcox.test(as.numeric(Monocyte)*100 ~ factor(severity.max), metadata_temp)
p4 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Monocyte)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 5)), colour = "black") + ggtitle("Monocyte")
wilcoxtest <- wilcox.test(as.numeric(T_NK)*100 ~ factor(severity.max), metadata_temp)
p5 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(T_NK)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 24)), colour = "black") + ggtitle("T/NK")
wilcoxtest <- wilcox.test(as.numeric(B)*100 ~ factor(severity.max), metadata_temp)
p6 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(B)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("B")
wilcoxtest <- wilcox.test(as.numeric(Plasmablast)*100 ~ factor(severity.max), metadata_temp)
p7 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Plasmablast)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 4)), colour = "black") + ggtitle("Plasmablast")
Figure S2E:
plot_grid(p1,p2,p3,p4,p5,p6,p7,ncol = 7)
We repeat this analysis on Day 0 samples only as there were no COVID-19-negative samples collected after Day 0:
metadata_temp <- metadata_filtered[metadata_filtered$Day == "D0",]
my.cols <- brewer.pal(3, "Set2")
wilcoxtest <- wilcox.test(as.numeric(Neutrophil_total)*100 ~ factor(COVID), metadata_temp)
p1 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Neutrophil_total)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("Total")
wilcoxtest <- wilcox.test(as.numeric(Mature_Neutrophil)*100 ~ factor(COVID), metadata_temp)
p2 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Mature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Mature")
wilcoxtest <- wilcox.test(as.numeric(Immature_Neutrophil)*100 ~ factor(COVID), metadata_temp)
p3 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Immature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Immature")
wilcoxtest <- wilcox.test(as.numeric(Monocyte)*100 ~ factor(COVID), metadata_temp)
p4 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Monocyte)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Monocyte")
wilcoxtest <- wilcox.test(as.numeric(T_NK)*100 ~ factor(COVID), metadata_temp)
p5 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(T_NK)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 4)), colour = "black") + ggtitle("T/NK")
wilcoxtest <- wilcox.test(as.numeric(B)*100 ~ factor(COVID), metadata_temp)
p6 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(B)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("B")
wilcoxtest <- wilcox.test(as.numeric(Plasmablast)*100 ~ factor(COVID), metadata_temp)
p7 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Plasmablast)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 15)), colour = "black") + ggtitle("Plasmablast")
Figure S3A:
plot_grid(p1,p2,p3,p4,p5,p6,p7,ncol = 7)
Next we want to understand determinants of sample purity. We can observe correlations between ANC and CIBERSORTx Total neutrophil fraction, as well as ALC and CIBERSORTx T/NK fraction.
metadata_temp <- metadata_filtered[metadata_filtered$COVID == "Positive",]
metadata_temp <- metadata_temp[complete.cases(metadata_temp$AMC.matched) & complete.cases(metadata_temp$ALC.matched),]
my.cols <- brewer.pal(3, "Set2")
p1 <- ggplot(metadata_temp, aes(x = factor(ANC.matched), y = (Neutrophil_total)*100)) + geom_boxplot(outlier.shape = NA, fill = "blue3") + geom_jitter(height = 0, alpha = 0.3) + theme_bw() + stat_compare_means() + ylab("CIBERSORTx Total Neutrophil (%)") + xlab("ANC Category") + scale_fill_manual(values = my.cols[c(1,2)]) + scale_y_continuous(limits = c(0,max(metadata_temp$Neutrophil_total)*100)) + coord_fixed(ratio = 1/14) + theme(legend.position = "none")
p2 <- ggplot(metadata_temp, aes(x = factor(ALC.matched), y = (T_NK)*100)) + geom_boxplot(outlier.shape = NA, fill = "mediumpurple4") + geom_jitter(height = 0, alpha = 0.3) + theme_bw() + stat_compare_means() + ylab("CIBERSORTx T/NK (%)") + xlab("ALC Category") + scale_fill_manual(values = my.cols[c(1,2)]) + scale_y_continuous(limits = c(0,max(metadata_temp$T_NK)*100)) + coord_fixed(ratio = 1/11) + theme(legend.position = "none")
p3 <- ggplot(metadata_temp, aes(x = factor(AMC.matched), y = (Monocyte)*100)) + geom_boxplot(outlier.shape = NA, fill = "skyblue") + geom_jitter(height = 0, alpha = 0.3) + theme_bw() + stat_compare_means() + ylab("CIBERSORTx Monocyte (%)") + xlab("AMC Category") + scale_fill_manual(values = my.cols[c(1,2)]) + scale_y_continuous(limits = c(0,max(metadata_temp$T_NK)*100)) + coord_fixed(ratio = 1/11) + theme(legend.position = "none")
Figure S3B:
plot_grid(p1,p2,p3,ncol=3)
We see that lower purity samples tend to have lower ANC and higher ALC. The converse is also true.
metadata_lo <- metadata_filtered[order(metadata_filtered$Neutrophil_total),]
metadata_lo <- metadata_lo[1:round(nrow(metadata_lo)*.25),]
comparisondf <- matrix(0L, nrow = 25, ncol = 3)
colnames(comparisondf) <- c("Var1","Var2","Freq")
comparisondf[,1] <- c(rep(5,5),rep(4,5),rep(3,5),rep(2,5),rep(1,5))
comparisondf[,2] <- c(rep(1:5,5))
for (i in 1:nrow(comparisondf)){
comparisondf[i,3] <- sum(metadata_lo$ANC.matched == comparisondf[i,1] & metadata_lo$ALC.matched == comparisondf[i,2], na.rm = TRUE)
}
theme_nogrid <- function (base_size = 12, base_family = "") {
theme_bw(base_size = base_size, base_family = base_family) %+replace%
theme(panel.grid = element_blank())
}
my.cols <- brewer.pal(7, "Greens")
p1 <- ggplot(as.data.frame(comparisondf), aes(x = factor(Var1), y = Freq, fill = factor(Var2))) + geom_bar(position="stack", stat="identity") + scale_fill_manual(values = (my.cols[3:7])) + xlab("ANC Category") + ylab("Count") + theme_bw() + coord_fixed(ratio = .1)
p1$labels$fill <- "ALC Quintile"
metadata_hi <- metadata_filtered[order(metadata_filtered$Neutrophil_total),]
metadata_hi <- metadata_hi[round(nrow(metadata_hi)*.75):nrow(metadata_hi),]
comparisondf <- matrix(0L, nrow = 25, ncol = 3)
colnames(comparisondf) <- c("Var1","Var2","Freq")
comparisondf[,1] <- c(rep(5,5),rep(4,5),rep(3,5),rep(2,5),rep(1,5))
comparisondf[,2] <- c(rep(1:5,5))
for (i in 1:nrow(comparisondf)){
comparisondf[i,3] <- sum(metadata_hi$ANC.matched == comparisondf[i,1] & metadata_hi$ALC.matched == comparisondf[i,2], na.rm = TRUE)
}
my.cols <- brewer.pal(7, "Greens")
p2 <- ggplot(as.data.frame(comparisondf), aes(x = factor(Var1), y = Freq, fill = factor(Var2))) + geom_bar(position="stack", stat="identity") + scale_fill_manual(values = (my.cols[3:7])) + xlab("ANC Category") + ylab("Count") + theme_bw() + coord_fixed(ratio = .1)
p2$labels$fill <- "ALC Quintile"
metadata_filtered$Purity <- -1*as.numeric(rownames(metadata_filtered) %in% rownames(metadata_lo)) + as.numeric(rownames(metadata_filtered) %in% rownames(metadata_hi))
metadata_filtered$Purity <- mapvalues(metadata_filtered$Purity, from = c(-1,0,1), to = c("lo","mid","hi"))
metadata_temp <- metadata_filtered[metadata_filtered$Purity %in% c("lo","hi"),]
fisher.test(table(metadata_temp$ANC.matched, metadata_temp$Purity))
##
## Fisher's Exact Test for Count Data
##
## data: table(metadata_temp$ANC.matched, metadata_temp$Purity)
## p-value < 2.2e-16
## alternative hypothesis: two.sided
fisher.test(table(metadata_temp$ALC.matched, metadata_temp$Purity))
##
## Fisher's Exact Test for Count Data
##
## data: table(metadata_temp$ALC.matched, metadata_temp$Purity)
## p-value = 0.01286
## alternative hypothesis: two.sided
Figure S3C:
p1
p2
Next, we view the CIBERSORTx estimated neutrophil fractions over time in COVID-19-positive patients:
my.cols <- brewer.pal(3, "Set2")
metadata_temp <- metadata_filtered[metadata_filtered$COVID == "Positive" & metadata_filtered$Day %in% c("D0","D3","D7"),]
kwtest <- kruskal.test(g = metadata_temp$Day, x = metadata_temp$Neutrophil_total)
p1 <- ggplot(metadata_temp, aes(x = factor(Day), y = as.numeric(Neutrophil_total)*100, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(kwtest$p.value, digits = 15)), colour = "black") + ggtitle("Total")
kwtest <- kruskal.test(g = metadata_temp$Day, x = metadata_temp$Mature_Neutrophil)
p2 <- ggplot(metadata_temp, aes(x = factor(Day), y = as.numeric(Mature_Neutrophil)*100, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(kwtest$p.value, digits = 37)), colour = "black") + ggtitle("Mature")
kwtest <- kruskal.test(g = metadata_temp$Day, x = metadata_temp$Immature_Neutrophil)
p3 <- ggplot(metadata_temp, aes(x = factor(Day), y = as.numeric(Immature_Neutrophil)*100, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(kwtest$p.value, digits = 34)), colour = "black") + ggtitle("Immature")
Figure 1D:
plot_grid(p1,p2,p3,ncol=3)
We also check if CIBERSORTx estimated neutrophil-related cell type fractions are associated with any clinical parameters.
my.cols <- brewer.pal(5,"Reds")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Creatinine.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p1 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Creatinine.matched), y = Neutrophil_total, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p2 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Creatinine.matched), y = Neutrophil_total, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p3 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Creatinine.matched), y = Neutrophil_total, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Greens")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CRP.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p4 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CRP.matched), y = Neutrophil_total, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p5 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CRP.matched), y = Neutrophil_total, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p6 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CRP.matched), y = Neutrophil_total, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Blues")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Ddimer.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p7 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Ddimer.matched), y = Neutrophil_total, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p8 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Ddimer.matched), y = Neutrophil_total, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p9 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Ddimer.matched), y = Neutrophil_total, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Purples")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$LDH.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p10 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(LDH.matched), y = Neutrophil_total, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p11 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(LDH.matched), y = Neutrophil_total, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p12 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(LDH.matched), y = Neutrophil_total, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
# my.cols <- brewer.pal(10,"Paired")
# metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$intubated),]
# metadata_temp <- metadata_temp[metadata_temp$COVID == "1",]
# p13 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(intubated), y = Neutrophil_total, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p14 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(intubated), y = Neutrophil_total, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p15 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(intubated), y = Neutrophil_total, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
my.cols <- brewer.pal(10,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CXR.infiltrates),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p16 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CXR.infiltrates), y = Neutrophil_total, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p17 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CXR.infiltrates), y = Neutrophil_total, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p18 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CXR.infiltrates), y = Neutrophil_total, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
my.cols <- brewer.pal(12,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Trop_72h),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p19 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Trop_72h), y = Neutrophil_total, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p20 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Trop_72h), y = Neutrophil_total, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p21 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Trop_72h), y = Neutrophil_total, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
Figure S4 (Total Neutrophil):
plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,ncol=3)
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
plot_grid(#p13,p14,p15,
p16,p17,p18,p19,p20,p21,ncol=3)
my.cols <- brewer.pal(5,"Reds")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Creatinine.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p1 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Creatinine.matched), y = Mature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p2 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Creatinine.matched), y = Mature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p3 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Creatinine.matched), y = Mature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Greens")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CRP.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p4 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CRP.matched), y = Mature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p5 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CRP.matched), y = Mature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p6 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CRP.matched), y = Mature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Blues")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Ddimer.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p7 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Ddimer.matched), y = Mature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p8 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Ddimer.matched), y = Mature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p9 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Ddimer.matched), y = Mature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Purples")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$LDH.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p10 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(LDH.matched), y = Mature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p11 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(LDH.matched), y = Mature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p12 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(LDH.matched), y = Mature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
# my.cols <- brewer.pal(10,"Paired")
# metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$intubated),]
# metadata_temp <- metadata_temp[metadata_temp$COVID == "1",]
# p13 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(intubated), y = Mature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p14 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(intubated), y = Mature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p15 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(intubated), y = Mature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
my.cols <- brewer.pal(10,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CXR.infiltrates),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p16 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CXR.infiltrates), y = Mature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p17 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CXR.infiltrates), y = Mature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p18 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CXR.infiltrates), y = Mature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
my.cols <- brewer.pal(12,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Trop_72h),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p19 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Trop_72h), y = Mature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p20 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Trop_72h), y = Mature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p21 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Trop_72h), y = Mature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
Figure S4 (Mature Neutrophil):
plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,ncol=3)
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
plot_grid(#p13,p14,p15,
p16,p17,p18,p19,p20,p21,ncol=3)
my.cols <- brewer.pal(5,"Reds")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Creatinine.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p1 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Creatinine.matched), y = Immature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p2 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Creatinine.matched), y = Immature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p3 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Creatinine.matched), y = Immature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Greens")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CRP.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p4 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CRP.matched), y = Immature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p5 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CRP.matched), y = Immature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p6 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CRP.matched), y = Immature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Blues")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Ddimer.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p7 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Ddimer.matched), y = Immature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p8 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Ddimer.matched), y = Immature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p9 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Ddimer.matched), y = Immature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
my.cols <- brewer.pal(5,"Purples")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$LDH.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p10 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(LDH.matched), y = Immature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p11 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(LDH.matched), y = Immature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p12 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(LDH.matched), y = Immature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
# my.cols <- brewer.pal(10,"Paired")
# metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$intubated),]
# metadata_temp <- metadata_temp[metadata_temp$COVID == "1",]
# p13 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(intubated), y = Immature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p14 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(intubated), y = Immature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p15 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(intubated), y = Immature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
my.cols <- brewer.pal(10,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CXR.infiltrates),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p16 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CXR.infiltrates), y = Immature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p17 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CXR.infiltrates), y = Immature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p18 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CXR.infiltrates), y = Immature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
my.cols <- brewer.pal(12,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Trop_72h),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p19 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Trop_72h), y = Immature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p20 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Trop_72h), y = Immature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p21 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Trop_72h), y = Immature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
Figure S4 (Immature Neutrophil):
plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,ncol=3)
## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
plot_grid(#p13,p14,p15,
p16,p17,p18,p19,p20,p21,ncol=3)
We next generate UMAPs for sample visualization in Figure 1E. Figure S5 features projections of many clinical parameters and quality control parameters - as is the case with other parameters in the study, not all can be reproduced here.
set.seed(10101)
#tpm.umap <- umap(t(logTPM_filtered))
#metadata_filtered$umap1 <- tpm.umap$layout[,1]
#metadata_filtered$umap2 <- tpm.umap$layout[,2]
metadata_filtered$severity.covid <- matrix(0L, nrow = nrow(metadata_filtered), ncol = 1)
for (i in 1:nrow(metadata_filtered)){
if (metadata_filtered$COVID[i] == "Positive"){
metadata_filtered$severity.covid[i] <- metadata_filtered$severity.max[i]
}
if (metadata_filtered$COVID[i] == "Negative"){
if (metadata_filtered$severity.max[i] == "H"){
metadata_filtered$severity.covid[i] <- "H"
}
else {
metadata_filtered$severity.covid[i] <- "non_COVID"
}
}
}
my.cols <- brewer.pal(3, "Set2")
p1 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2)) + theme_bw() + geom_point(aes(colour = factor(COVID)),size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + scale_color_manual(values = c(my.cols[1], my.cols[2]))
p1$labels$colour <- "COVID"
my.cols.2 <- brewer.pal(3, "RdBu")
p2 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2)) + theme_bw() + geom_point(aes(colour = factor(severity.covid)),size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + scale_color_manual(values = c("grey", my.cols[1], my.cols.2[3], my.cols[2]))
p2$labels$colour <- "Severity"
a <- "Mature_Neutrophil"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered[,colnames(metadata_filtered) == a])*100,max(metadata_filtered[,colnames(metadata_filtered) == a])*100))
p3 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2, colour = (as.numeric(Mature_Neutrophil)*100))) + theme_bw() + geom_point(size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + sc
p3$labels$colour <- "Mature Neu"
a <- "Immature_Neutrophil"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered[,colnames(metadata_filtered) == a])*100,max(metadata_filtered[,colnames(metadata_filtered) == a])*100))
p4 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2, colour = (as.numeric(Immature_Neutrophil)*100))) + theme_bw() + geom_point(size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + sc
p4$labels$colour <- "Immature Neu"
Figure 1E:
plot_grid(p1,p2,p3,p4,ncol=2)
We perform differential expression to test for genes and pathways associated with COVID-19 infection. We use the helper script Neutrophil_DESeq2.R to import the function. As we will see, if we do not make a correction for plasmablast contamination, we will have many highly expressed immunoglobulin genes in the COVID+ side.
source(paste0(prefix,"Neutrophil_DESeq2.R"))
rownames(metadata_filtered) <- metadata_filtered$Public.Sample.ID
DESeq2_list <- Neutrophil_DESeq2(counts = Count_filtered, mdata = metadata_filtered, day = "D0")
dds <- DESeqDataSetFromMatrix(countData = DESeq2_list$Count_select, colData = DESeq2_list$coldata, design = ~ COVID)
dds <- DESeq(dds)
res <- as.data.frame(results(dds, name="COVID_Positive_vs_Negative"))
filenam <- "Day0_COVIDpos_vs_COVIDneg_uncorrected"
temp <- genepc[which(genepc$Gene.stable.ID %in% rownames(res)),]
res$symbol <- matrix(0L, nrow = nrow(res))
for (i in 1:nrow(res)){
if (rownames(res)[i] %in% temp$Gene.stable.ID){
res$symbol[i] <- temp$Gene.name[which(rownames(res)[i] == temp$Gene.stable.ID)]
} else {
res$symbol[i] <- rownames(res)[i]
}
}
res$rank <- sign(res$log2FoldChange)*(-1)*log10(res$pvalue)
res <- res[complete.cases(res),]
res_sig <- res[res$padj < 0.05,]
#write.table(res,paste0(prefix,"DESeq2/",filenam,".txt"),sep = "\t")
resordered <- res[order(res$rank),]
log2fc <- as.numeric(resordered$log2FoldChange)
log10p <- as.numeric(-1*log10(resordered$pvalue))
padj <- resordered$padj
rank <- resordered$rank
symbol <- resordered$symbol
combo <- cbind(log2fc,log10p,padj,symbol,rank)
colnames(combo) <- c("log2fc","log10p","padj","symbol","rank")
rownames(combo) <- rownames(resordered)
combo <- as.data.frame(combo)
combo$rank <- as.numeric(combo$rank)
combo$log10p <- as.numeric(combo$log10p)
combo$log2fc <- as.numeric(combo$log2fc)
combo$significance <- as.numeric(combo$padj < 0.05)
combo$significance <- as.factor(combo$significance)
IGgenes <- c("IGKC","IGKV4-1","IGKV5-2","IGKV6-21","IGKV3D-20","IGKV3D-11","IGLV4-69","IGLV8-61","IGLV4-60","IGLV6-57","IGLV10-54","IGLV1-51","IGLV1-47","IGLV7-46","IGLV5-45","IGLV1-44","IGLV7-43","IGLV1-40","IGLV1-36","IGLV3-27","IGLV3-25","IGLV2-23","IGLV3-21","IGLV3-19","IGLV2-18","IGLV3-16","IGLV2-14","IGLV2-11","IGLV3-10","IGLV3-9","IGLV3-1","IGLC1","IGLC2","IGLC3","IGLC7","IGHA2","IGHG4","IGHG2","IGHA1","IGHG1","IGHG3","IGHD","IGHM","IGHV6-1","IGHV1-2","IGHV1-3","IGHV2-5","IGHV3-7","IGHV3-11","IGHV3-13","IGHV3-15","IGHV1-18","IGHV3-20","IGHV3-21","IGHV3-23","IGHV1-24","IGHV2-26","IGHV3-33","IGHV4-34","IGHV4-39","IGHV1-46","IGHV3-48","IGHV3-49","IGHV5-51","IGHV3-53","IGHV1-58","IGHV4-61","IGHV3-66","IGHV1-69","IGHV2-70D","IGHV3-73","IGLV9-49","IGHV3-64","IGKV3D-15","IGHV4-59","IGHV3-74","IGHV3-72","IGHV4-31","IGHV3-43","IGKV2D-30","IGKV1-6","IGKV3-20","IGKV1D-33","IGKV1-17","IGKV1-8","IGKV1-16","IGKV1D-16","IGKV2-24","IGKV3-11","IGKV1-9","IGKV1-33","IGKV1-39","IGKV2D-28","IGKV1D-17","IGKV2-30","IGKV2D-29","IGKV1-12","IGKV1-5","IGKV2-28","IGKV3-15","IGKV1-27","IGKV2D-40","IGKV1D-39","IGLVI-70","IGKV2-29","IGLL5","IGHV3-30","IGHV2-70","IGKV1D-13","IGHV4-4","IGLV2-8","IGHV1-69D","IGHV7-4-1","IGHV3-64D","IGHV5-10-1")
combo$color <- 0
combo$color[combo$symbol %in% IGgenes] <- 1
combo$labels <- 0
combo$labels[combo$symbol %in% c("MZB1","JCHAIN")] <- 1
combo$labels <- as.factor(combo$labels)
options(ggrepel.max.overlaps = Inf)
my.cols <- brewer.pal(3,"Set2")
plot1 <- ggplot(combo, aes(x = log2fc, y = log10p)) + geom_point(data = subset(combo, color == 1), colour = my.cols[2]) + geom_point(data = subset(combo, color == 0), colour = "grey") + geom_text_repel(data = subset(combo, labels == 1), aes(label = as.character(symbol))) + theme_bw() + ylab("-Log10(p-value)") + xlab("Log2(FC)") + annotate("text", x=1.3, y=0, label= "COVID+", colour = my.cols[2]) + annotate("text", x=-1.2, y=0, label= "COVID-", colour = my.cols[1]) + coord_fixed(ratio = .13) + theme(panel.grid = element_blank())
Figure S6A:
plot1
Thus we will use CIBERSORTx estimated cell type fractions as well as an additional immunoglobulin score to regress non-neutrophil contamination out of the data.
source(paste0(prefix,"Pathway_scoring.R"))
IG.score <- Pathway_scoring(IGgenes)
metadata_filtered$IG.score <- (IG.score)
We check the correlations with this score with MZB1, a plasmablast marker, and MS4A1, a B cell marker.
goi <- "MZB1"
metadata_filtered$GOI <- t(logTPM_filtered[genepc$Gene.stable.ID[which(genepc$Gene.name == goi)],])
metadata_filtered$GOI <- as.numeric(metadata_filtered$GOI)
summry <- lm(GOI ~ IG.score, metadata_filtered)
p1 <- ggplot(metadata_filtered, aes(x = as.numeric(IG.score), y = GOI)) + geom_point() + theme_bw() + geom_smooth(method = "lm") + xlab("Immunoglobulin Score") + ylab("MZB1 Log2(TPM+1)") + annotate("text", x=-3, y=8.75, label= paste0("R2 =",round(summary(summry)$r.squared,3)), colour = "black") + annotate("text", x=-3, y=8, label= paste0("p =",round(summary(summry)$coefficients[,4],308)), colour = "black")
goi <- "MS4A1"
metadata_filtered$GOI <- t(logTPM_filtered[genepc$Gene.stable.ID[which(genepc$Gene.name == goi)],])
metadata_filtered$GOI <- as.numeric(metadata_filtered$GOI)
summry <- lm(GOI ~ IG.score, metadata_filtered)
p2 <- ggplot(metadata_filtered, aes(x = as.numeric(IG.score), y = GOI)) + geom_point() + theme_bw() + geom_smooth(method = "lm") + xlab("Immunoglobulin Score") + ylab("MS4A1 Log2(TPM+1)") + annotate("text", x=-3, y=5, label= paste0("R2 =",round(summary(summry)$r.squared,3)), colour = "black") + annotate("text", x=-3, y=4.5, label= paste0("p =",round(summary(summry)$coefficients[,4],18)), colour = "black")
Figure S6B:
p1
p2
View the immunoglobulin score on a UMAP:
a = "IG.score"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered[,colnames(metadata_filtered) == a]),max(metadata_filtered[,colnames(metadata_filtered) == a])))
p1 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2, colour = (as.numeric(IG.score)))) + theme_bw() + geom_point(size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + sc
p1$labels$colour <- "IG Score"
Figure S6C:
p1
We check the IG score across COVID status, Day, and SeverityMax.
my.cols <- brewer.pal(3,"Set2")
p1 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = IG.score, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.15) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .5) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + stat_compare_means()
Figure S6D:
p1
my.cols <- brewer.pal(5,"Set2")
p2 <- ggplot(metadata_filtered, aes(x = factor(Day), y = IG.score, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.15) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .5) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3],my.cols[4],my.cols[5])) + stat_compare_means()
Figure S6E:
p2
my.cols <- brewer.pal(3,"RdBu")
p3 <- ggplot(metadata_filtered[metadata_filtered$severity.max %in% c("severe","non-severe"),], aes(x = factor(severity.max), y = IG.score, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.15) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .5) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + stat_compare_means()
Figure S6F:
p3
Now we can perform the COVID+ vs. COVID- differential expression analysis.
DESeq2_list <- Neutrophil_DESeq2(counts = Count_filtered, mdata = metadata_filtered, day = "D0")
dds <- DESeqDataSetFromMatrix(countData = DESeq2_list$Count_select, colData = DESeq2_list$coldata, design = ~ Neutrophil_total + T_NK_factor + Monocyte_factor + IG_factor + Plasmablast_factor + COVID)
dds <- DESeq(dds)
res <- as.data.frame(results(dds, name="COVID_Positive_vs_Negative"))
filenam <- "Day0_COVIDpos_vs_COVIDneg_correct-NeuCont+TNK+Monocyte+Plasmablast+IG"
temp <- genepc[which(genepc$Gene.stable.ID %in% rownames(res)),]
res$symbol <- matrix(0L, nrow = nrow(res))
for (i in 1:nrow(res)){
if (rownames(res)[i] %in% temp$Gene.stable.ID){
res$symbol[i] <- temp$Gene.name[which(rownames(res)[i] == temp$Gene.stable.ID)]
} else {
res$symbol[i] <- rownames(res)[i]
}
}
res$rank <- sign(res$log2FoldChange)*(-1)*log10(res$pvalue)
res <- res[complete.cases(res),]
res_sig <- res[res$padj < 0.05,]
#write.table(res,paste0(prefix,"DESeq2/",filenam,".txt"),sep = "\t")
resordered <- res[order(res$rank),]
log2fc <- as.numeric(resordered$log2FoldChange)
log10p <- as.numeric(-1*log10(resordered$pvalue))
pvalue <- as.numeric(resordered$pvalue)
padj <- resordered$padj
rank <- resordered$rank
symbol <- resordered$symbol
combo <- cbind(log2fc,log10p,padj,pvalue,symbol,rank)
colnames(combo) <- c("log2fc","log10p","padj","pvalue","symbol","rank")
rownames(combo) <- rownames(resordered)
combo <- as.data.frame(combo)
combo$rank <- as.numeric(combo$rank)
combo$log10p <- as.numeric(combo$log10p)
combo$log2fc <- as.numeric(combo$log2fc)
combo$pvalue <- as.numeric(combo$pvalue)
combo$significance <- as.numeric(combo$padj < 0.05)
combo$significance <- as.factor(combo$significance)
combo$color <- 0
combo$color[combo$pvalue < 1e-4 & combo$log2fc > 0.5] <- 1
combo$color[combo$pvalue < 1e-4 & combo$log2fc < -0.5] <- -1
combo$labels <- 0
combo$labels[rownames(combo) %in% rownames(resordered)[1:20]] <- 1
combo$labels[rownames(combo) %in% rev(rownames(resordered))[1:20]] <- 1
combo$labels <- as.factor(combo$labels)
options(ggrepel.max.overlaps = Inf)
my.cols <- brewer.pal(3,"Set2")
plot1 <- ggplot(combo, aes(x = log2fc, y = log10p)) + geom_point(data = subset(combo, color == 0), colour = "grey") + geom_point(data = subset(combo, color == 1), colour = my.cols[2]) + geom_point(data = subset(combo, color == -1), colour = my.cols[1]) + geom_text_repel(data = subset(combo, labels == 1), aes(label = as.character(symbol))) + theme_bw() + ylab("-Log10(p-value)") + xlab("Log2(FC)") + annotate("text", x=1.3, y=0, label= "COVID+", colour = my.cols[2]) + annotate("text", x=-1.2, y=0, label= "COVID-", colour = my.cols[1]) + coord_fixed(ratio = .21) + theme(panel.grid = element_blank())
Figure 1F:
plot1
Next we perform Gene Set Enrichment Analysis.
gmt.file <- gmtPathways(paste0(prefix,"all_gene_sets.gmt"))
res <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 14)
ranking <- res[,"rank"]
names(ranking) <- res$symbol
set.seed(15001)
fgseaRes <- fgsea(pathways = gmt.file,
stats = ranking,
minSize=25,
maxSize=1000,
eps = 0)
#write.table(fgseaRes[,1:7], file = paste0(prefix,"GSEA_",filenam,".txt"), sep = "\t")
The running enrichment score data points are generated using the following code.
getEnrichmentDataframe <- function (pathway, stats, gseaParam = 1, ticksSize = 0.2) {
rnk <- rank(-stats)
ord <- order(rnk)
statsAdj <- stats[ord]
statsAdj <- sign(statsAdj) * (abs(statsAdj)^gseaParam)
statsAdj <- statsAdj/max(abs(statsAdj))
pathway <- unname(as.vector(na.omit(match(pathway, names(statsAdj)))))
pathway <- sort(pathway)
gseaRes <- calcGseaStat(statsAdj, selectedStats = pathway, returnAllExtremes = TRUE)
bottoms <- gseaRes$bottoms
tops <- gseaRes$tops
combo <- as.data.frame(cbind(tops,bottoms))
combo$average <- matrix(0L, nrow = nrow(combo), ncol = 1)
for (p in 1:nrow(combo)){
combo$average[p] <- (combo$tops[p]+combo$bottoms[p])/2
}
combo <- cbind(combo,pathway)
return(combo)
}
signalingpathways <- c("HALLMARK_INTERFERON_GAMMA_RESPONSE","HALLMARK_INTERFERON_ALPHA_RESPONSE","GO_CYTOKINE_MEDIATED_SIGNALING_PATHWAY","HALLMARK_TNFA_SIGNALING_VIA_NFKB","GO_INTERLEUKIN_1_BETA_PRODUCTION","GO_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY","HALLMARK_IL6_JAK_STAT3_SIGNALING","GO_INTERLEUKIN_8_PRODUCTION","GO_TUMOR_NECROSIS_FACTOR_SUPERFAMILY_CYTOKINE_PRODUCTION","GO_INTERLEUKIN_6_PRODUCTION","GO_INSULIN_LIKE_GROWTH_FACTOR_RECEPTOR_SIGNALING_PATHWAY","GO_RESPONSE_TO_EPIDERMAL_GROWTH_FACTOR")
cellularpathways <- c("GO_DEFENSE_RESPONSE_TO_VIRUS","GO_INNATE_IMMUNE_RESPONSE","GO_REGULATION_OF_VIRAL_GENOME_REPLICATION","GO_OXIDATIVE_PHOSPHORYLATION","GO_HYDROGEN_PEROXIDE_CATABOLIC_PROCESS","GO_REGULATION_OF_CARBOHYDRATE_BIOSYNTHETIC_PROCESS","GO_RIBOSOME_ASSEMBLY","GO_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE")
dataframe <- getEnrichmentDataframe(gmt.file[[signalingpathways[1]]], ranking)
signalingdf <- cbind(dataframe$pathway,dataframe$average,rep(signalingpathways[1],nrow(dataframe)))
colnames(signalingdf) <- c("rank","enrichment","pathway")
for (i in 2:length(signalingpathways)){
dataframe <- getEnrichmentDataframe(gmt.file[[signalingpathways[i]]],
ranking)
temp <- cbind(dataframe$pathway,dataframe$average,rep(signalingpathways[i],nrow(dataframe)))
colnames(temp) <- c("rank","enrichment","pathway")
signalingdf <- rbind(signalingdf,temp)
}
dataframe <- getEnrichmentDataframe(gmt.file[[cellularpathways[1]]], ranking)
cellulardf <- cbind(dataframe$pathway,dataframe$average,rep(cellularpathways[1],nrow(dataframe)))
colnames(cellulardf) <- c("rank","enrichment","pathway")
for (i in 2:length(cellularpathways)){
dataframe <- getEnrichmentDataframe(gmt.file[[cellularpathways[i]]],
ranking)
temp <- cbind(dataframe$pathway,dataframe$average,rep(cellularpathways[i],nrow(dataframe)))
colnames(temp) <- c("rank","enrichment","pathway")
cellulardf <- rbind(cellulardf,temp)
}
my.cols <- brewer.pal(12,"Paired")
p1 <- ggplot(as.data.frame(signalingdf), aes(x = as.numeric(rank), y = as.numeric(enrichment), colour = pathway)) + geom_point() + theme_bw() + scale_colour_manual(values = rev(my.cols)) + coord_fixed(ratio = 10000) + theme(panel.grid.minor = element_blank(), legend.text=element_text(size=6)) + xlab("Rank in Gene List") + ylab("Running Enrichment Score") + ggtitle("GSEA: Signaling Pathways")
my.cols <- brewer.pal(8,"Dark2")
p2 <- ggplot(as.data.frame(cellulardf), aes(x = as.numeric(rank), y = as.numeric(enrichment), colour = pathway)) + geom_point() + theme_bw() + scale_colour_manual(values = rev(my.cols)) + coord_fixed(ratio = 10000) + theme(panel.grid.minor = element_blank(), legend.text=element_text(size=6)) + xlab("Rank in Gene List") + ylab("Running Enrichment Score") + ggtitle("GSEA: Cellular Processes")
Figure 1G:
p1
Figure 1H:
p2
sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] fgsea_1.18.0 umap_0.2.7.0
## [3] ggpubr_0.4.0 heatmap3_1.1.9
## [5] SeuratObject_4.0.2 Seurat_4.0.4
## [7] cowplot_1.1.1 ggmosaic_0.3.3
## [9] corrplot_0.90 openxlsx_4.2.4
## [11] DESeq2_1.32.0 SummarizedExperiment_1.22.0
## [13] Biobase_2.52.0 MatrixGenerics_1.4.3
## [15] matrixStats_0.60.1 GenomicRanges_1.44.0
## [17] GenomeInfoDb_1.28.2 IRanges_2.26.0
## [19] S4Vectors_0.30.0 BiocGenerics_0.38.0
## [21] dplyr_1.0.7 plyr_1.8.6
## [23] RColorBrewer_1.1-2 ggrepel_0.9.1
## [25] ggplot2_3.3.5 knitr_1.33
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.20 tidyselect_1.1.1
## [4] RSQLite_2.2.8 AnnotationDbi_1.54.1 htmlwidgets_1.5.3
## [7] grid_4.1.1 BiocParallel_1.26.2 Rtsne_0.15
## [10] munsell_0.5.0 codetools_0.2-18 ica_1.0-2
## [13] future_1.22.1 miniUI_0.1.1.1 withr_2.4.2
## [16] colorspace_2.0-2 highr_0.9 ROCR_1.0-11
## [19] ggsignif_0.6.2 tensor_1.5 listenv_0.8.0
## [22] labeling_0.4.2 GenomeInfoDbData_1.2.6 polyclip_1.10-0
## [25] bit64_4.0.5 farver_2.1.0 parallelly_1.27.0
## [28] vctrs_0.3.8 generics_0.1.0 xfun_0.25
## [31] fastcluster_1.2.3 R6_2.5.1 locfit_1.5-9.4
## [34] bitops_1.0-7 spatstat.utils_2.2-0 cachem_1.0.6
## [37] DelayedArray_0.18.0 assertthat_0.2.1 promises_1.2.0.1
## [40] scales_1.1.1 gtable_0.3.0 globals_0.14.0
## [43] goftest_1.2-2 rlang_0.4.11 genefilter_1.74.0
## [46] splines_4.1.1 rstatix_0.7.0 lazyeval_0.2.2
## [49] spatstat.geom_2.2-2 broom_0.7.9 yaml_2.2.1
## [52] reshape2_1.4.4 abind_1.4-5 backports_1.2.1
## [55] httpuv_1.6.2 tools_4.1.1 ellipsis_0.3.2
## [58] spatstat.core_2.3-0 ggridges_0.5.3 Rcpp_1.0.7
## [61] zlibbioc_1.38.0 purrr_0.3.4 RCurl_1.98-1.4
## [64] rpart_4.1-15 openssl_1.4.5 deldir_0.2-10
## [67] pbapply_1.4-3 zoo_1.8-9 haven_2.4.3
## [70] cluster_2.1.2 magrittr_2.0.1 data.table_1.14.0
## [73] RSpectra_0.16-0 scattermore_0.7 lmtest_0.9-38
## [76] RANN_2.6.1 fitdistrplus_1.1-5 hms_1.1.0
## [79] patchwork_1.1.1 mime_0.11 evaluate_0.14
## [82] xtable_1.8-4 XML_3.99-0.7 rio_0.5.27
## [85] readxl_1.3.1 gridExtra_2.3 compiler_4.1.1
## [88] tibble_3.1.4 KernSmooth_2.23-20 crayon_1.4.1
## [91] htmltools_0.5.2 mgcv_1.8-36 later_1.3.0
## [94] tidyr_1.1.3 geneplotter_1.70.0 DBI_1.1.1
## [97] MASS_7.3-54 Matrix_1.3-4 car_3.0-11
## [100] igraph_1.2.6 forcats_0.5.1 pkgconfig_2.0.3
## [103] foreign_0.8-81 plotly_4.9.4.1 spatstat.sparse_2.0-0
## [106] annotate_1.70.0 XVector_0.32.0 stringr_1.4.0
## [109] digest_0.6.27 sctransform_0.3.2 RcppAnnoy_0.0.19
## [112] spatstat.data_2.1-0 Biostrings_2.60.2 rmarkdown_2.10
## [115] cellranger_1.1.0 leiden_0.3.9 fastmatch_1.1-3
## [118] uwot_0.1.10 curl_4.3.2 shiny_1.6.0
## [121] lifecycle_1.0.0 nlme_3.1-152 jsonlite_1.7.2
## [124] carData_3.0-4 viridisLite_0.4.0 askpass_1.1
## [127] fansi_0.5.0 pillar_1.6.2 lattice_0.20-44
## [130] KEGGREST_1.32.0 fastmap_1.1.0 httr_1.4.2
## [133] survival_3.2-13 glue_1.4.2 zip_2.2.0
## [136] png_0.1-7 bit_4.0.4 stringi_1.7.4
## [139] blob_1.2.2 memoise_2.0.0 irlba_2.3.3
## [142] future.apply_1.8.1